Optimal control for non-Markovian open quantum systems 
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An efficient optimal control theory based on the Krotov method is introduced for a non-Markovian 
open quantum system with a time-nonlocal master equation in which the control parameter and the 
bath correlation function are correlated. This optimal control method is developed via a quantum 
dissipation formulation that transforms the time-nonlocal master equation to a set of coupled linear 
time-local equations of motion in an extended auxiliary Liouville space. As an illustration, the 
optimal control method is applied to find the control sequences for high-fidelity Z-gates and identity- 
gates of a qubit embedded in a non-Markovian bath. Z-gates and identity-gates with errors less 
than 10 -5 for a wide range of bath decoherence parameters can be achieved for the non-Markovian 
open qubit system with control over only the a z term. The control-dissipation correlation, and the 
memory effect of the bath are crucial in achieving the high-fidelity gates. 

PACS numbers: 03.65.Yz, 02.30.Yy, 03.67.Pp, 03.67.-a 



I. INTRODUCTION 

Quantum optimal control theory (QOCT) lh7i] is a 
powerful tool that provides a variational framework for 
calculating the optimal shaped pulse to maximize a de- 
sired physical objective (or minimize a physical cost func- 
tion). It has been applied to various open quantum sys- 
tems or models to obtain control sequences for quantum 
gate operations |8rU3|. Compared to the dynamical- 
decoupling-based method [14| in which a succession of 
short and strong pulses designed to suppress decoher- 
ence is applied to the system, QOCT is a continuous 
dynamical modulation with many degrees of freedom for 
selecting arbitrary shapes, durations and strengths for 
time-dependent control, and thus allows significant re- 
duction of the applied control energy and the correspond- 
ing quantum gate error. The authors of Ref. [8| investi- 
gated the optimal control of a qubit coupled to a two- 
level system that is exposed to a Markovian heat bath. 
Although this may mimic the reduced non-Markovian 
dynamics of the qubit, it is by no means a model of a 
qubit coupled directly to a non-Markovian environment. 
The authors of Refs. |9J-[ll| investigated optimal quan- 
tum gate operations in the presence of non-Markovian 
environments. However, to combine QOCT with a non- 
Markovian master equation involving time-ordered inte- 
gration of the nonunitary (dissipation) terms for noncorn- 
muting system and control operators, and for a nonlocal- 
in-time memory kernel, the numerical treatment is rather 
mathematically involved and computationally demand- 
ing. All of the QOCT approaches mentioned above [g- 
Il3| | for open quantum systems employed gradient-based 
[I Q algorithms for optimization. 

A somewhat different QOCT approach from the stan- 
dard gradient optimization methods is the Krotov it- 
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erative method @, [a, 0, UM ■ The Krotov method has 
several appealing advantages @, H, 0] over the gradient 
methods: (a) monotonic increase of the objective with 
iteration number, (b) no requirement for a line search, 
and (c) macrosteps at each iteration. A version of the 
Krotov optimization method has been used recently in 
Ref. [16| to deal with the non-Markovian optimal control 
problem of a quantum Brownian motion model with an 
exact stochastic equation of motion (master equation). 
We note however that only a few non-Markovian open 
quantum system models are exactly solvable |16l - t23| (the 
quantum Brownian model is one of them), and the ex- 
act master equations of these exactly solvable models are 
known to be in a time-local (time-convolutionless) form 
with time-dependent decoherence or decay rates, without 
involving the time-ordering problem of non-commuting 
operators. |16h23| . With time- local equations of mo- 
tion, the Krotov method could be directly employed to 
deal with optimal-control problems. Although it is com- 
mendable to derive an exact master equation, not too 
many problems can be exactly worked out in this way. 
For example, no exact master equation can be derived 
for the non-Markovian qubit-environment (spin-boson) 
model studied as an illustration for quantum gate opera- 
tions in this paper. It is thus important that an efficient 
QOCT approach based on the Krotov method for per- 
turbative treatment of general (not limited to just some 
certain classes of) non-Markovian open quantum systems 
should be developed. The perturbative non-Markovian 
master equation under only the Born approximation (or 
the weak system-bath coupling approximation) Il7j ] is 
in the form of a time-ordered non-commuting integro- 
differential equation [e.g., see Eqs. (P) and ©], and at 
first sight it is not at all clear how to effectively combine 
this kind of master equation with the Krotov optimiza- 
tion method. The study presented in this paper provides 
just such an efficient QOCT approach based on the Kro- 
tov method to deal with time-nonlocal non-Markovian 



open quantum systems. Our approach transforms the 
time-ordered non-commuting intcgro-diffcrcntial master 
equation into a set of time-local coupled differential equa- 
tions with the small price of introducing auxiliary den- 
sity matrices in an extended auxiliary Liouville space. As 
a result, incorporation of the resultant time-local equa- 
tions with the Krotov optimization method becomes ef- 
fective. We then apply the developed Krotov method to 
the problem of finding the optimal quantum gate control 
sequence for a qubit in a non-Markovian environment 
(bath). Our results illustrate that the control parameter 
can be engineered to efficiently counteract and suppress 
the environment effect for non-Markovian open systems 
with long bath correlation times (long memory effects). 
We also find that high-fidelity quantum gates with error 
smaller than 10 -5 can be achieved at long gate operation 
times for a wide range of bath decoherence parameters. 
This is in contrast to the cases in the literature |Sl— Tlljj . 
where the non-Markovian systems were mainly studied 
in a parameter regime very close to Markovian systems, 
and thus no significant reduction of the quantum gate 
errors was observed. 



II. NON-MARKOVIAN MASTER EQUATION 
AND QUANTUM OPTIMAL CONTROL THEORY 

In realistic experiments, one may have only limited 
control over the system Hamiltonian, and the maximum 
control parameter strength that can be realized is also 
restricted. Let us consider a total system with Hamilto- 
nian H = Hs(i) + Hj + Hb- Here the system Hamil- 
tonian Hs(t) = —he(t)<T z /2 — hfla x /2, describes a qubit 
with a time-dependent control parameter e(t) and a fixed 
tunneling frequency CI, i.e., having control over only the 
a z term. The bath Hamiltonian is Hb = J"^ huj q b q 'b q , 

where b q (b q ) is the creation (annihilation) operator of 
the bath mode q with frequency uj q . The interaction 
Hamiltonian Hi between the system and the bath with- 
out making the rotating-wave approximation is of the 
form Hj — cr x ^2 c q (b q + b q '), where c q is the coupling 
constant of bath mode q to the qubit system. Follow- 
ing the standard perturbation theory, we obtain (see the 
Appendix n for the derivation) under only the Born ap- 
proximation the time-convolution master equation for the 
reduced system density matrix as |24H27| 



p(t) = C s (t)p(t) + {C x K(t) + [£*£(<)]+}. (1) 

Here the superoperators Cs(t) and C x are defined via 
their actions on an arbitrary operator A, respectively, as 
C s {t)A = ±[H s (t),A] and C X A = &[a x ,A}. The non- 
Hermitian (dissipation) operator can be written as 



K{t) 



ih 



dt'C(t-t')Us{t,t')a xP {t'), 



(2) 



ing the time-ordering operator, and the bath correla- 
tion function (CF) is C(t — t 1 ) = J Q duJ(uj) cos[u(t — 

t')] coth[nw/(2k B T)] - i J °° du J(u) sin[w(t - t')] with T 
being the temperature. Note that Eq. ^ contains the 
bath CF C(t,t') and the time-ordered system propaga- 
tor superoperator lis{t, t') which involves the control pa- 
rameter e(t) through Hg{t) in Cs(t). Thus the control 
parameter and bath-induced nonunitary (dissipation) ef- 
fect are correlated. This paves the way for manipulat- 
ing the control sequence to counteract the effect of the 
bath on the system dynamics. This coherent control of 
non-Markovian open quantum systems is in contrast to 
various Markovian approaches of engineering reservoirs 
[28| and incoherent controls by directly manipulating the 
environments |29| . 

In the framework of the QOCT, one would like to max- 
imize the quality (fidelity) value of some target at time 
tf . Suppose the desired state-independent unitary quan- 
tum gate operation is denoted as O and the target is 
to perform a state-independent quantum gate operation 
as close as possible to O. We choose the trace distance 
between the desired target superoperator O and the ac- 
tual nonunitary propagator superoperator X(tf) at the 
final operation time tf to characterize the gate error, i.e., 
(error) = Tr{[£>- X(t f )} 2 }/N, where N is the dimension 
value of O. As minimizing the trace distance is similar 
to maximizing the real part of the trace fidelity [8(, we 
choose T = Re[Tr{C>t^( i/ )}]/7V as a quality (fidelity) 
measure of how well X(tf) approaches the target O in 
the QOCT framework. In realistic control problems, it is 
desirable that the optimal control sequence can provide 
highest quality (fidelity) with minimum energy consump- 
tion. Thus we introduce an objective function of the form 



J = T- 



dt , X(t'){e(t')-e (t')f 



(3) 



where the unitary qubit system propagator superop- 
erator U s (t,t') = T + exp{/*dr£ s (r)} with T+ be- 



where X(t) is a positive function that can be adjusted and 
chosen empirically, e(t) is the control parameter and £o(t) 
is a reference control value that can be properly chosen 
[2|, y] . Then the task of quantum control is to maximize 
the objective ([3]) under the constraint of the equation of 
motion of X{t) obtained by replacing p(t) with X(t) in 
Eqs. (HJ and ©. 

The time-dependent control parameter s(t) enters into 
the exponent of the time-ordered system propagator 
Us(t,t') 7 and thus appears inside the memory kernel or 
dissipation operator Eq. ([2]). As a result, Eq. (TTJ is a non- 
local time-ordered integro-differential equation in which 
values of the control parameter at all earlier times come 
into play, and is difficult to incorporate within the frame- 
work of QOCT. Thus an approach that retains the mer- 
its of the Krotov optimization method and removes the 
time-ordering and nonlocal problems in the equation of 
motion is very much desired. 



III. MASTER EQUATION AND OPTIMAL 
CONTROL IN EXTENDED LIOUVILLE SPACE 

One important observation to deal with the time- 
nonlocal non-Markovian quantum master equation is to 
express the bath CF's in a multi-exponential form [2J- 
E3, C(t-t') = £\ C J (0)e^(*-*'), where C,-(0) and lj are 
complex constants and can be found by numerical meth- 
ods. Then Eq. ([2]) can be written as JC(t) = 5D ■ /Cj- (i) , 

where £,-(*) - ^ p dt l C 3 {<d)e^^ t - t '^U s {t-t')(j x p(t'). Al- 
though A^j (£) is still a time- nonlocal and time-ordered in- 
tegration for non-commuting operators, if one takes the 
time derivative of Kj(t), one obtains 

t 3 {t) = (l/»fi)C,(0)cr x p(t) + [£<?(£) + 7i]/C,-(t), (4) 

with initial condition /Cj(0) = 0. The same process can 
be done for the Hermitian conjugate lO(t) = J2j %-j ($)■ 
Equation (JT]) combined with Eq. (J4|) and its Hermitian 
conjugate form a set of coupled linear equations of mo- 
tion that can be written as pit) = K(t)p\t), in terms of 
pit) = {pit) , JCj , K,p ; j = 1,2,3,...} in an extended aux- 
iliary Liouville space |24l - [27j . Obviously, the above equa- 
tions are time-local and have no time-ordering and mem- 
ory kernel integration problems. This yields a simple, fast 
and stable iterative scheme to incorporate with the Kro- 
tov method. The formal solution of p(t) can be written as 
pit) = G(t)f>i0), where the associated propagator super- 
operator can be shown to satisfy [dQit)/dt\ = A(t)(/(£) 
with 5(0) = Ijsf. Here I/j is the identity operator in 
the extended Liouville space and J\f is the dimension of 
Git). The real part of the trace fidelity for the propa- 
gator $(tf) is T = Re[Ti{QWit f )}}/J\f, where Q is the 
target operator O in the extended Liouville space. The 
goal of quantum optimal control here is to reach a de- 
sired target Q with maximum objective function J (or 
fidelity T) in a certain time t f. The optimal algorithm 
following the Krotov method [15j is summarized as fol- 
lows @, 0,0- (i) Guess an initial control sequence e it). 
(ii) Use the equations of motion to find the forward prop- 
agator Qkit) with the initial condition G(0) = 7v (fc = 
for the first iteration), (iii) Find an auxiliary backward 
propagator Bu{t) with the condition B(tf) = Q) . (iv) 
Propagate Gk+i(t) again forward in time, and update 
the control parameter iteratively with the rule £k+i(t) — 
£ kit) + 2^jRe[Tr{4(t)^gg fc+1 (t)}]. (v) Repeat steps 
(iii) and (iv) until a preset fidelity (error) is reached or 
until a given number of iterations has been performed. 



IV. RESULTS AND DISCUSSIONS 

To find the optimal control sequence for a state- 
independent single-qubit gates in a non-Markovian bath, 
one needs to know the bath spectral density to calculate 
the bath CF. We can deal with any form of the bath spec- 
tral density, but for simplicity, we take an Ohmic spectral 
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FIG. 1. (Color online) (a) Unitary Z-gate error versus op- 
eration time tf. (b) Bloch vector trajectory of an optimal 
unitary Z-gate operation with the initial polarization condi- 
tion P(0) = [1,0,0]. The inset in (a) is the optimal control 
pulse sequence for any tf > 0.3Q,. 
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FIG. 2. (color online) (a) Z-gate error versus time tf for 
different values of a and T at u c = 20$!. The stopping crite- 
rion of the error threshold is set to 10 -6 or when the number 
of iterations exceeds 3000. (b) Corresponding bath CF's for 
a = 0.01. We set fcs = 1 and h — 1 in all the figures presented 
here and below. 



density in the form of J(w) = auje^^^", where a is a di- 
mensionless damping constant, and w c is the bath cutoff 
frequency. The bath CF can thus be calculated and can 
then be expanded directly in a multi-exponential form. 
The values of Cj(0) and 7* of the exponentials are ob- 
tained numerically with the requirement that the error 
between the actual CF and the approximated CF is cho- 
sen to be less than or equal to 1CU 7 . Only three or four 
exponential terms in the expansion are required to ex- 
press the bath CF with such high accuracy as compared 
to a sum of more than 48 exponentials needed to express 
the same bath CF at a low temperature of T = 0.2J1 
through the spectral density parametrization [2J|. We 
note here that both the parametrization of Ref. [24| . 
which is highly efficient at moderate to high tempera- 
tures, and the direct bath CF decomposition of our ap- 
proach are well-suited for the parametrization of highly 
structured spectral densities, leading to long and oscil- 
latory bath CF's. However, our approach has a great 
computational advantage over the commonly used ap- 
proach of spectral density parametrization |24| - l27j at low 
temperatures. 

The ideal Z-gate performance as a function of the gate 
operation time tf in the absence of the bath is given in 
Fig. Q] with the stopping criterion of error set to 10 ~ 8 . 
For all the calculations performed in this paper, we re- 




FIG. 3. (Color online) (a) Z-gate error versus time tf for 
the same bath parameters and stopping criterion as in Fig. [2] 
except uj c — CI. (b) Corresponding bath CF's for a = 0.1. 
The lower three curves are plotted in a smaller vertical-axis 
range in the inset. 




FIG. 4. (Color online) (a) Z-gate error versus cutoff ui c for 
different values of a and T with tf = SI" 1 , (b) Optimal Z- 
gate control sequence for a = 0.01, T — Cl, tf — CI" 1 and 



CI. 



strict the maximum control parameter e(t) < 3057. Ex- 
cellent Z-gate performance can be achieved for pulse time 
tf > 0.3/57. The corresponding optimal pulse sequence 
is shown in the inset of Fig. HJa). In fact, a unitary Z- 
gate for pulse time tf > 0.3/57 with error limited only by 
machine precision can be achieved. Furthermore, if one 
does not impose any restriction on the maximum control 
parameter strength, an optimal perfect Z-g&te can be 
achieved for any finite period of time tf. This is in con- 
trast to the control pulse strategy for the implementation 
of a (unitary) Z-gate in [8J , where the gate operation time 
tf > 7r/S7, corresponding to the static Cla x inducing at 
least a full loop around the x axis. One can see from the 
Bloch polarization vector evolution in Fig. [IJb) that our 
optimal control pulse strategy does not require a full loop 
around the x axis as compared with the state evolution 
in Fig. 2 (left panel) of [8j, and thus can achieve a Z-gate 
in a much shorter operation time limited primarily by 
the control parameter strength. 

Figures [2] and [3] show the errors of the Z-gate versus 
operation times and the corresponding bath CF's for dif- 
ferent values of the bath dimensionless damping constant 
a and temperature T at cutoff frequencies lo c = 2057 and 
lu c = 57, respectively. Note that the imaginary part of 
the bath CF does not depend on temperature T, so only 
one curve of the imaginary part is plotted in Figs. [5Jb) 
and[3Jb). One can see that for uj c = 2057, the gate er- 
rors increase with the operation times tf for all the dif- 



ferent bath parameter cases studied in Fig. Ufa). The 
errors for the case of smaller a = 0.01 are about 10~ 3 at 
the beginning and increase to 10 _1 at longer operation 
times. As expected, the gate performance deteriorates 
as the value of a increases. The value of a plays a more 
important role than that of temperature in determining 
the amount of the error for the u> c = 2057 case. The 
bath CF at time t = for the high-temperature case of 
ksT = 10fi-57 is bigger than that for the low temperature 
case of fcgT = M7 by a factor of only 1.5 as shown in 
Fig. &b). However, the bath CF (effective decay rate) 
is directly proportional to a. Therefore, the gate errors 
for the case u> c = 2057 shown in Fig. [2ja) depend mainly 
on the value of a. In all the previous investigations of 
quantum gate performance for both Markovian and non- 
Markovian open quantum systems [8l4ll| , the gate errors 
for a value of a = 0.01 are about 10~ 3 or worse. This is 
similar to what we found for the u> c — 2057 case. How- 
ever, one of our main results, significantly different from 
previous investigations, is that for small w c , it is possi- 
ble to achieve high-fidelity Z-gates with errors smaller 
than 10~ 5 at large times, as shown in Fig.[3Ja). This can 
be understood from the control-dissipation correlation, 
Eq. ([2]) , and the memory effect of the non-Markovian en- 
vironment. The bath CF's in Figs. Ufa) and^a) decay to 



zero on a time scale of ts 



called the bath correla- 



tion time. For uj c = 2057, the bath is very close to Marko- 
vian, and the dissipator ([2]) or decay rate approaches a 
constant value in a very short time of ts ~ 0.0557" 1 
and thus cannot be significantly suppressed by the ex- 
ternal control sequence. On the other hand, for uj c = 57, 
the longer bath correlation time of ts ~ 57 _1 allows the 
optimal control sequence to have enough action time to 
counteract and suppress the contribution from the bath 
in the weak system-bath coupling case, and thus reduces 
the gate error considerably. For a — 0.01 and temper- 
ature ksT < hCl, the gate error can be kept about or 
smaller than 10 -6 . The temperature T plays a similar 
role to the coupling strength when ksT > M7, as the 
amplitude ratio of the bath CF between the fcgT = 10M7 
and ksT < Kl cases is about 10. Even though the gate 
errors for larger a and higher T increase with tf, the 
errors for all of the cases studied except for the case of 
a = 0.1 and T = 1057 in Fig.[^a) are smaller than 10" 5 , 
smaller than the error threshold of 10 -4 (10 -3 ) [30] re- 
quired for fault-tolerant quantum computation. We also 
perform calculations for an optimal identity gate and the 
behaviors of the gate errors versus operation times tf are 
similar to that of the Z-gate in Fig. [21a). Achieving a 
high-fidelity identity gate at long times implies having 
the capability for arbitrary state preservation, i.e., stor- 
ing arbitrary state robustly against the bath. The gate 
errors are expected to be much lower if one has indepen- 
dent control over both the a z and a x terms in the qubit 
Hamiltonian and if there is no restriction on the control 
parameter strengths. 

Figure H)Ja) shows the Z-gate error versus bath cutoff 
frequency w c at operation time tf = 57 _1 for different 



values of a and T. One can see that the gate error de- 
pends strongly on the bath cutoff frequency. The error 
increases as w c becomes bigger. For the weak coupling 
and low temperature cases (a = 0.01, T < Q), it is pos- 
sible to reduce the error to below 10~ 5 for ui c < 2.50. 
These indicate that the bath correlation time ts ~ w^T 1 
or memory effect plays an important role in determining 
the gate error. Figure 0|b) shows the optimal control se- 
quence for uj c — CI, a = 0.01, ksT = Ml and tf = £l~ l . 
The optimal control sequence that suppresses the deco- 
herence induced by the bath is totally different from that 
of the ideal unitary case in the inset of Fig. [Ha). 



noise spectroscopy techniques [3lM35| to determine the 
bath spectral densities experienced by the qubits and 
then applying QOCT to find explicitly control sequences 
for quantum gate operations for the non-Markovian open 
qubit system will yield fast quantum gates with low in- 
vested energy and high fidelity . By virtue of its general- 
ity and efficiency, our Krotov based QOCT method will 
find useful applications in many different branches of the 
sciences. Recent experiments on engineering external en- 
vironments [40J, simulating open quantum systems [41|, 
and observing non-Markovian dynamics |42l |43j could fa- 
cilitate the experimental realization of the QOCT in non- 
Markovian open quantum systems in the near future. 



V. CONCLUSION 

To conclude, a universal QOCT based on the Krotov 
method for a time-nonlocal non-Markovian open quan- 
tum system has been introduced and applied to ob- 
tain control sequences and gate errors for Z and iden- 
tity gates. Our study has yielded several computa- 
tional and conceptual innovations: (a) The QOCT that 
combines the Krotov method and an extended Liouville 
space quantum dissipation formulation that transforms 
the nonlocal-in-time master equation to a set of cou- 
pled linear local-in-time equations of motion is intro- 
duced to deal with non-Markovian open quantum sys- 
tems, (b) The direct numerical decomposition of the 
bath CF into multi-exponential form has a great com- 
putational advantage over the commonly used spectral 
density parametrization approach [2J-|27j. (c) The con- 
structed QOCT, which retains the merits of the Kro- 
tov method, is extremely efficient in dealing with the 
time-nonlocal non-Markovian equation of motion. Com- 
pared to the calculations performed on a 40-node SUN 
Linux cluster via the gradient-based approach to tackle 
the nonlocal kernel directly Q, the calculations using 
our approach for a similar problem can be performed 
on a typical laptop PC with ease, thus opening the way 
for investigating two-qubit and many-qubit problems in 
non-Markovian environments, (d) Our study of opti- 
mal control reveals the strong dependence of the gate 
errors on the bath correlation time and exploits this non- 
Markovian memory effect for high-fidelity quantum gate 
implementation and arbitrary state preservation in an 
open quantum system. The presented QOCT has been 
shown to be a powerful tool, capable of facilitating imple- 
mentations of various quantum information tasks against 
decoherence. The required information is knowledge of 
the bath or noise spectral density which is experimentally 
accessible by, for examp le, dyn amical decoupling noise 
spectroscopy techniques [3lM35| . We note here that not 
only our proposed method of QOCT but also dynamical 
decoupling and other strategies for fighting decoherence 
require knowledge of the spectral distribution of the noise 
(bath spectral density) in order to improve the strategies 
and design effective and/or optimized control sequences 
[12l Il3l . 1361-39] . Thus using the dynamical decoupling 
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Appendix: Derivation of the time-nonlocal master 
equation 

Here we provide the derivation of the time-nonlocal 
non-Markovian master equation (fTJ) and © in the text. 
Following the assumption of factorized initial system- 
bath state pr(0) — Pt(0) — p(0) <8> Ro, the standard 
perturbative time-convolution (under only the Born ap- 
proximation) master equation in the interaction picture 
takes the form of 



p(t) = —rTr B [H I (t),p(0)®R o ] 
n 



-^Tr B / df[fr I (t),[i? / (t') >/ 5(t / )®-Ro]IA.l) 

Here p(i) is the system density matrix in the interaction 
picture, Rq = exp(— i/s/A;BT)/TrB[exp(— iJa/fc^T)] is 
the initial thermal reservoir density operator at temper- 
ature T, and the system-bath interaction Hamiltonian in 
the interaction picture in our spin-boson model can be 
written as 



Hi(t) = a x {t)B{t), (A.2) 

where a x {t) = Ul(t)a x lls(t) with U s (t) = 
T + ei Jo Hs ( t ) dt being the system evolution opera- 
tor and T + being the time-ordering operator, and 

&Je 1 ""*. Substituting Eq. (fAT2|) 



B(t)=EgC q b q e~ ,u ' 

into Eq. (|A.1[) and then tracing over the reservoir (bath 



-1"q 



or environment) degrees of freedom, one obtains 



n 2 



K 2 



P(t)= "Is I ^ B [d x (t)B(t),[a x {t')B(t'), p{t') <g> R ]]dt' 

dt' {[a x (t)a x (t')p(t') 

-a x (t')p(t')a x (t)]C(t-t') 
+\p(t')a x (t')a x (t) 

-a x {t)p{t')a x {t')]C{t' -*)} , (A.3) 

where the relation Ti B [Hj(t)R ] — has been used to 
eliminate the first-order term in Eq. (IA.1I) . The bath 
correlation function is defined as 

C{t-t')=Tr B [B(t)B(t')R } 

dwJ(o;){[n(a;) + l]e^(*- t ') 

+n(u J )e lu ' {t - t ' ) \ 
i 
duiJ(ui) cos[ui(t — t')] coth 



hcu 



2k B T 



du)J(ui) s'm[u)(t — t')], 



(A.4) 



where we have used the definition of B(t) below Eq. (|A.2|) 
and the definition of the spectral density J(w) = 
J2 q \ c q\ 2 ^{°° ~ w <?) to obtain Eq. (jA.4p , and n{ui) is the 



canonical ensemble average occupation number of the 
bath. Rotating back to the Schrodinger picture, the 
nonlocal-in-time non-Markovian master equation for the 
system density matrix p(t) takes the form 



m = 

where 



l[H s (t),p(t)]-'-{[a x ,)C(t)]-P(t),a x ]}, 



(A.5) 



-U s (t) 



[ dt' 
Jo 



C(t-t')a x (t')p(t') 



ul(t), 

(A, 



and pit) = U s (t)p(t)Us(t). We can further rewrite 
Eqs. (|A.5[) and (|A.6[) in a superoperator form as 



P (t) = c s (t)p(t) + {c x ic(t) + [c x ic(t)]*}, 



(A.7) 



£(*) = -*/ dt'C(t-t')U s (t,t')a x p(t'). (A.8) 



and 



Here the superoperators Cs{t) and C x are defined via 
their actions on an arbitrary operator A, respectively, as 
C s {t)A = -(i/H)[H s (t),A] and C X A = -(i/H)[a x ,A], 
and the system propagator superoperator Us(t,t') = 
T + exp{/ t ! drCsir)}. Equations (|A.7|) and (|A.8|) are just 
Eqs. (Ti]) and © in the text. 
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